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Technical  Objectives 

The  research  supported  by  ONR  grant  number  NOOO 14- 1 3-1-0665  was  conducted  in  two  phases,  as  the 
initial  stages  of  this  effort  involved  closing  significant  research  on  the  subsurface  signature  of  wake- 
radiated  internal  waves.  The  primary  objectives  of  each  phase  are  listed  below. 

STAGE-l,  “Subsurface  Signature  of  the  Internal  Wave  Field  Radiated  by  Submerged  High  Reynolds 
Number  Stratified  Wakes”: 

•  Characterize  the  low-wind-condition  subsurface  signature  of  a  single  high-frequency  internal  wave 
beam(s)  in  a  uniform  or  vertically  variable  stratification,  the  latter  typical  of  littora]  ocean  conditions. 

•  Correlate  qualitatively  and  quantitatively,  over  a  range  of  Reynolds  and  Froude  numbers,  the  patterns 
of  subsurface  motion  observed  above  a  stratified  turbulent  wake  with  the  flow  patterns  inside  the 
wake  core  to  infer  the  state  (“age”)  of  wake  turbulence. 

STAGE-2,  “Reynolds  number  scaling  of  Stratified  Turbulent  Wakes”: 

•  Identify  the  Reynolds  number  scaling  for  Non-Equilibrium  (NEQ)  regime  duration,  diagnostics  of 
secondary  turbulence,  vertical  momentum  flux  and  radiated  internal  wave  properties  in  a  stratified 
turbulent  wake. 

•  Develop  Reynolds/Froude  number  parameterizations  of  associated  vertical  transport  coefficients  for 
use  in  RANS  and  self-similarity-based  models  ;  assess  the  accuracy  of  the  latter  by  comparing  with 
LES  results. 

Note  that,  prior  to  its  originally  planned  termination  for  12/31/2016,  on  account  of  changes  in  the  ONR 
logistical  systems,  this  project  was  replaced  in  late  Spring  2015  by  award  N00014-1 5-1-25 13,  titled 
“High  Reynolds  Number  Stratified  Turbulent  Wakes:  Internal  Wave  Energetics,  Seif-Similarity  and 
Subgrid-Scafe  Modeling”,  The  objectives  of  this  new  award  essentially  overlap  with  and  extend  those 
listed  in  STAGE-2  above. 


Results 


The  results  obtained  in  pursuit  of  the  above  objectives  may  be  grouped  in  the  following  categories,  on 
each  of  which  we  further  elaborate  below; 

(i)  Lagrangian  flows  in  the  reflection  of  an  internal  wave  beam  off  a  free-slip  surface. 

(ii)  Surface/subsurface  signature  of  the  wake-radiated  internal  wave  field. 

(iii)  High  Reynolds  number  stratified  wakes. 

(iv)  Collaborative  DNS/LES  of  high  Reynolds  number  stratified  wakes  (Frontier-award  support). 


(i)  Lagrangian  flows  in  tfie  reflection  of  an  internal  wave  beam  off  a  free-slip 
surface. 

A  significant  part  of  the  previous  grant  involved  the  quantification  of  mean  flows  in  the  reflection  region 
of  an  internal  wave  beam  (IWB)  incident  on  a  free-slip  surface,  for  a  linear  stratification  either  extending 
to  the  top  or  with  a  mixed  layer  overlaying  it  (see  schematic  in  Figure  I).  Results  on  mean  flow  scaling 
and  parameterization  and  comparison  with  theoretical  findings  of  a  group  in  MIT  were  published 
previously  by  Zhou  and  Diamessis  (Phys.  Fluids  2013)  through  previous  ONR  support.  However,  a 
perplexing  question  was  repeatedly  posed  to  the  PI  and  PhD  student  at  conferences:  Are  these  mean 
flows  Eulerian  or  Lagrangian  in  nature  ?  Any  wave-driven  flow  has  associated  with  it  a  Stokes  drift,  i.e. 
the  velocity  linked  to  the  residual  displacement  of  a  particle  after  it  has  completed  one  full  cycle  of  its 
orbit  over  one  wave  period.  The  Stokes  drift  is  very  likely  to  counterbalance  the  Eulerian  mean  flow,  and 
the  true  Lagrangian  mean  may  be  dramatically  different  than  what  our  previous  study  proposed. 


Figure  1.  Schematic  of  the  reflection  of  an  internal  wave  beam  (TWB)  off  a  free-slip  surface.  The 
beam  is  incident  at  an  angle  0  to  the  vertical.  The  shaded  region  is  the  location  where  nonlinear 
dynamics  and,  therefore,  Lagrangian  mean  drift  are  most  potent. 

An  extensive  existing  database  of  19  2-D  simulations  of  free-slip  reflecting  IWBs  in  linear  stratifications 
was  to  quantify  the  Lagrangian  flows  in  the  reflection  region.  The  parameters  typically  varied  were  the 
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IWB  amplitude,  inclination  9  to  the  vertical  and  beam  compactness  (wavelengths  contained  in  the  beam 
envelope). 

Straightforward  theoretical  analysis  shows  that  an  0(A)  IWB  beam  (here  ^  is  a  measure  of  the  beam 
amplitude,  typically  given  by  its  steepness,  i.e.,  max.  isopycnal  displacement  normalized  by  horizontal 
wavelength)  generates  an  0{A“)  mean  flow  and  an  0(A^)  first  harmonic  in  the  reflection  zone.  Moreover, 
one  can  show  motions  of  0(A^)  can  produce  a  0(A^^^)  Stokes  drift.  In  particular,  analytical  estimation  of 
the  Stokes  drift  in  the  imiscid  theoretical  surrogate  of  our  modeled  IWB  indicates  that  it  exactly  cancels 
our  the  Eulerian  mean  flow  for  all  beams  considered.  In  other  words,  at  least  on  a  theoretical  level,  no 
mean  transport  of  surface  particulates  should  be  expected  in  the  (sub)surface  zone. 

Nevertheless,  the  theoretical  beam  neglects  both  viscous  effects  but  also  higher-order  effects.  Hence,  we 
performed  Lagrangian  particle  tracking  on  all  of  the  above  19  datasets.  The  particle  tracking  algorithm  is 
based  on  spectral  interpolation  in  space  in  the  horizontal  (Fourier)  and  vertical  (multidomain-based 
Legendre/Lagrange).  Additionally,  it  employs  an  OiAt"^)  Adams-Bashforth-Moulton  scheme  in  time.  The 
high  spatiotemporal  accuracy  is  paramount  to  produce  reliable  results  over  the  long-time  integrations 
associated  with  the  monitoring  of  each  particle.  Parallel  implementation  of  the  algorithm  has  been 
enabled  through  use  of  par-for  loops  in  Matlab. 


Figure  2,  Typical  type~A  particle  orbits  from  from  TOl  {A  =  3,23%  and  45*^).  Trajectories  of  a 
two-dimensional  array  of  16^16  —  256  particles  are  plotted  over  five  wave  periods.  Left  panel 
shows  the  full  field  of  view  spanned  by  the  original  array  of  particles  ;  right  panel  shows  a  zoomed- 
in  view  within  the  dashed-line  box  drawn  in  the  left  panel.  The  drifts  in  the  mean  orbital  position 
are  so  small  that  orbits  from  five  periods  effectively  overlap  into  one  closed  trajectory, 
overlap  onto  one,  for  each  initial  position. 


Within  the  reflection  zone,  over  a  region  of  dimension  A-  (X^  and  X.  are  the  wavelengths  in  the 
horizontal  and  vertical  direction,  respectively)  square  array  of  16><16  particles  is  inserted  over  18  equally 
spaced  intervals  within  one  wave  period.  This  total  of  of  4,068  particles  is  tracked  over  the  course  of  five 
IWB  periods  generating  the  desired  Lagrangian  statistics. 
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We  have  classified  our  simulations  into  three  broader  categories,  based  on  the  observed  particle  orbits 
within  the  reflection  region  (Figures  2  to  4).  Type-A  trajectories,  typical  of  lower-amplitude  IWBs 
oriented  closer  to  the  vertical  are  fairly  regular  and  closed,  i.e.,  no  net  drill  of  the  particles  is  observed. 
As  the  beam  orientation  is  moved  away  further  from  the  vertical,  at  an  inclination  angle  0=63° ,  a  weak 
drift  is  observed  in  the  particle  position  of  the  order  of  0.0  U,.  Such  particle  orbits  are  regarded  as  Type- 
B.  At  the  same  value  of  6  and  with  only  a  twofold  increase  in  IWB  amplitude,  the  drift  becomes  stronger 
and  particle  orbits,  regarded  as  Type-C,  are  subject  to  significant  distortions. 


■0.2 


-0.4 


-0.0 
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-0.6  *■ 
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Figure  3:  Typical  type-B  particle  orbits  from  from  TOl  (A  =  1.27%  and  0  =  63°).  A  net  drift  can 
now  be  observed  for  each  initial  position,  while  the  orbital  shape  remains  repeatable  from  one 
wave  cycle  to  the  next. 

Depth-averaged  values  of  the  mean  and  rms  IWB-mean  Lagrangian  drifts  (not  shown  here)  show  that  in 
the  mean,  the  net  drift  is  not  that  powerful  because  of  variations  in  its  sign  in  the  vertical  direction. 
However,  the  rms  value  of  the  Lagrangian  drift  is  significant,  reaching  as  high  as  10%  of  the  IWB  phase 
speed.  Figure  5  (left  panel)  shows  the  maximum  rms  value  of  the  Lagrangian  (sum  of  Euterian  and 
Stokes)  drift  velocity  Ut  as  a  function  of  .4tan0,  a  parameter  which  successfully  scaled  the  Eulerian  mean 
flows  investigated  in  our  previous  studies.  Unlike  the  theoretical  estimate,  C4  is  not  zero  but  scales  aA^ , 
although  it  is  much  weaker  than  its  Eulerian  contribution  (Figure  5)  with  which  it  follows  the  same 
power  law  behavior.  Moreover,  when  normalized  by  the  horizontal  phase  speed,  q,  it  can  get  reach 
values  as  high  as  10%. 

The  occurrence  of  strong  rms  values  of  the  Lagrangian  drift  vs.  significantly  lower  mean  values  of  it, 
suggests  that  although  there  might  not  be  a  collective  horizontal  drift  of  particles  across  the  water 
column,  significant  particle  dispersion  will  indeed  take  place.  This  dispersion  may  be  quantified  in  the 
form  of  a  horizontal  dispersion  coefficient,  Kx,  defined  by  o^=2KxAt  where  a  is  the  variance  of  particle 
displacements  over  a  characteristic  time  At,  taken  to  be  one  or  multiple  wave  period(s).  Straightforward 
arguments  show  that  scales  as  A*  as  is  indeed  confirmed  by  our  data  for  At  equal  to  3  wave  periods 
(Figure  5  right  panel). 


4 


Ti2:  Initial  phase  no.  1 


Figure  4:  Typical  type-C  particle  orbits  from  from  TOl  (^4  =  2,65%  and  0  =  63^),  Particle  orbits  are 
shown  for  two  wave  cycles  only;  Trajectories  during  the  first  cycle  are  drawn  in  blue,  and  the 
immediately  following  cycle  in  red.  Particles  are  downsampled  in  the  vertical  for  clarity*  The 
orbital  shape  of  some  particles  varies  from  the  previous  wave  cycle  to  the  next,  and  a  great  amount 
of  net  drift  in  the  mean  orbital  position  can  be  observed. 
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Figure  5;  Left  panel:  Quadratic  scaling  ViJcx^{Ai2inOy  of  Lagrangian  mean  drift  velocity  scale  £4* 
Dash  line  is  the  quadratic  fit  for  all  data  points  of  against  At&nB;  dash-dot  line  is  the 
theoretical  scaling  of  Eulerian  mean  flow  <u>/ca  =  8;r^(Atantf)^  computed  in  previous  ONR-funded 
work.  Right  panel;  Quartic  scaling  -'{AtauOyaf  dispersion  coefficient  within  a 

reflecting  IWB,  Dash  line  represents  power  law  (Atan0)^ 

The  question  that  naturally  arises  is  how  one  might  relate  these  results  to  the  surface  reflection  of  wake- 
radiated  internal  waves,  under  at  least  the  common  assumption  that  an  IWB  is  a  reliable  first-order 
surrogate  of  such  waves  ?  In  past  work,  we  have  found  that  the  internal  waves  generated  by  stratified 
sphere  wakes  have  horizontal  wavelengths  in  the  range  [jD,3£)]  for  4  <  F  r<  64.  The  wave  orientation 
angle  increases  up  to  0=55^  for  higher  Re  wakes,  with  these  waves  potentially  assuming  an  even  more 
horizontal  inclination  at  Re  values  not  accessible  to  parallel  LES  nowadays.  For  these  same  waves  we 


5 


have  also  found  that  c,  «  0.0 1 U,  where  U  is  the  sphere  tow  speed.  Extrapolating  the  above  to  operational 
conditions,  where  Z>=10m  and  U~  lOm/s,  this  suggests  a  Lagrangian  drift  velocity  of  {4*  0.001{/=  1 
cm/s  .  Given  that  internal  wave  periods  are  typically  a  factor  of  two  longer  than  the  ambient  buoyancy 
period  Af',  with  typically  0.01<A^<0.001  Hz  one  might  expect  drifts  of  5  to  50m  over  five  buoyancy 
periods.  Our  dispersion  coefficient  estimates  suggest  that  these  drifts  will  be  accompanied  by  average 
particle  displacements  of  cr  =  =  IQT^ XJ"  =~O(10m). 

Our  findings  on  Lagrangian  drifts  and  dispersion  caused  by  a  reflecting  internal  wave  beam  may  be 
found  in  the  article  by  Zhou  and  Diamessis  (Phys.  Fluids  2015). 

(ii)  Surface/subsurface  signature  of  the  wake-radiated  internal  wave  field. 


An  extensive  analysis  of  six  LES  datasets  in  very  deep  and  wide  domains  was  performed  to  explore  the 
signature  of  wake-radiated  internal  waves  on  an  rigid  lid  surface  within  an  idealized  linear  stratification. 
The  datasets  span  parameter  values  of  internal  Froude  number,  Fr=4,  16  and  64,  and  two  values  of 
/?e=5’<10^and  10^ .  The  problem  geometry  is  shown  in  Figure  6,  which  shows  a  schematic  of  an  Oyz 
transect.  The  stratification  used  is  always  linear,  as  more  complex  stratifications  (e.g.  pycnocline)  and 
background  currents  are  issues  that  will  considered  in  future  studies. 


Free-slip  rigid  lid 


Wave  absorbing  layer 


Figure  6:  Schematic  of  a  span-depth  cross-section  of  the  computational  domain  used  in  the  LES  of 
the  (sub)surface  signature  of  wake-radiated  IWs.  The  wake  core  is  located  at  the  same  depth  of  9D 
for  all  runs,  whereas  the  domain  is  sufficiently  wide  to  allow  for  the  surface  reflection  of  the  wake- 
emitted  waves  which  eventually  dissipated  in  the  absorbing  layers. 


A  fundamental  pillar  in  the  analysis  of  the  surface-reflecting  wake-emitted  internal  wave  field  are  1-D 
and  2-D  wavelets.  1-D  Morlet  wavelets  are  used  to  compute  the  dominant  internal  wave  frequency  as  a 
function  of  spanwise  offset  and  time.  2-D  Morlet  (non-directionally  biased)  wavelets  are  used  to  identify 
the  horizontal  internal  wave  lengths  at  the  surface.  Details  of  the  computation  procedure  are  reported  in 
the  thesis  of  the  recently  graduated  PhD  student  supported  by  this  grant. 
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2-D  wavelet  analysis  of  the  surface-reflecting  internal  wave  fields  exploiting  the  homogeneity  of  the 
streamwise  direction,  generates  a  number  of  data  points  in  a  3-D  volume  in  the  space  (y/D,  M) 
where  >/Z)  and  Xf/D  are  the  non-dimensional  spanwise  offset  and  horizontal  wavelength,  respectively. 
This  volume  of  data  may  be  collapsed  along  the  y  direction,  as  Ah  is  not  found  to  exhibit  a  strong 
transverse  dependence  {unlike  the  actual  frequency  which  is  typical  of  a  highly  dispersive  wave  field).  A 
joint  p.d.f*  in  the  Nt)  may  then  be  constructed,  as  shown  in  Figure  7»  The  peak  values  of  Joint 
p.dT/s  correspond  to  the  energetic  internal  waves  with  the  most  frequently  occurring  wavelength 
impacting  the  free  surface  and  straining  it  accordingly.  At  a  given  Re^  the  energetic  packets  arrive  at  the 
surface  at  a  later  time  as  Fr  is  reduced.  For  a  given  Fr,  the  arrival  time  is  shifted  in  time  towards  a  later 
M.  Finally,  the  most  frequently  occurring  wavelengths  appear  to  lie  along  a  universal  line  with  the 
location  of  the  centers  of  the  joint  p.d.f.  contours  varying  as  a  function  of  {Re.Fr), 


Figure  7:  Joint  probability'  density^  contours  of  most  energetic  wave  packets  in  the 
(Xj/DyNt)  sample  space  at  (a)  Re  =  S  and  (b)  Re  —  respectively  for 

Fr  e{4, 16,  64}. 


At  a  given  A’r,  one  can  average  along  A//A  for  values  of  probability  density  exceeding  a  certain 
threshold.  Figure  8  shows  the  evolution  of  the  resulting  conditionally-averaged  mean  horizontal 

wavelength  D  across  all  six  cases  considered.  The  evolution  of  follows  a  universal  trend,  for  all 
simulations  performed,  described  by  a  {Ntf^  power  law.  Such  a  power  law  is  consistent  with  the  findings 
of  previous  experimental  studies  and  motivates  one  to  revert  to  linear  theory. 

The  definition  of  the  vertical  group  velocity,  according  to  linear  theory,  ultimately  leads  to  a  relationship 
between  Xi^D  and  the  origin,  zo,  of  the  wake-radiated  internal  waves  and  their  angle  of  propagation, 

_  z  —  zq  27r  1 

D  D  sinflcos^  ^  iV(f to) 

Here  z=0  corresponds  to  the  wake  centerline.  If  one  assumes  that  t»tQ,  i*e.,  significant  time  has  elapsed 
since  the  waves  were  emitted  and,  when  impacting  the  free  surface,  they  are  in  the  'Tar-field”  w/r  to  the 
wake  core,  Xf/D  may  be  approximated  as: 
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z  —  zq  2y  1 
D  s.va9 cos^ 9  Nt' 
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Figure  $:  Time  series  of  conditionally  averaged  wavelengths  I D  observed  at  the  sea 
surface.  Best  fit  to  the  linear  propagation  model  is  shown  as  dash  line. 

Notice  how  a  (M)’’  power  law  is  recovered  from  purely  linear  theory.  This  finding  suggests  that  the  key 
process  driving  the  dynamics  of  the  internal  waves  upon  arrival  at  the  free  surface  is  linear  dispersion.  In 
other  words,  at  least  for  the  Re  and  Fr  considered  here,  the  far*field  behavior  of  the  wake-emitted  internal 
waves  is  linear ;  any  non-linearities  are  restricted  to  the  near-field  region,  immediately  around  the  wake 
edge.  Any  non-linear  effects  at  the  surface  originate  from  the  interaction  of  an  incident  packet  with  its 
reflecting  counterpart  and  not  any  interactions  among  upward  propagating  wave  packets. 

Can  any  prediction  be  made  on  the  depth  of  the  wave-generating  wake  based  on  surface  information, 
namely  measured  values  of  Ai/D  ?  To  this  end,  a  schematic  is  given  in  Figure  9.  Consider  the  surface 
depth  of  z=9£)  that  holds  for  all  the  simulations  considered  here.  The  surface  needs  to  be  placed 
sufficiently  far  from  the  wave-generating  turbulence  for  the  above  equation  to  hold,  i.e.,  the  waves  need 
to  propagate  over  a  sufficiently  long  distance  to  actually  disperse.  Now,  considering  a  typical  value  of 
^=45“  and  using  a  least-squares  fit  power  law  for  the  data  in  Figure  8,  one  finds  that  Z(/D  =  0.62  ±0.11  . 
The  “virtual  origin”  of  the  waves  should  not  be  confused  with  the  actual  wake  edge  and  is,  instead,  a 
calibration  parameter  of  our  linear  model,  much  like  a  virtual  origin  may  be  identified  for  a  wake  or  jet 
by  relying  on  the  corresponding  self-similarity  solution. 
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Figure  9;  Schematic  of  key  concepts  associated  with  our  linear  wave  model  for  wake-radiated 
internal  wave  propagation  as  evaluated  at  the  sea  surface  (z  =  Internal  waves  (dash  arrows)  of 
various  wavelengths  >1// are  observed  at  the  surface  as  if  they  had  been  emitted  from  the  virtual 
origin  Zo  at  4  «  L  The  linear  model  relates  the  surface  observed  to  the  submerged  turbulence 
characterized  by  the  dimensionless  timeM  since  the  passage  of  the  wake-generating  body^ 

Nevertheless,  returning  to  the  original  question  above,  inferring  M  solely  from  Xh/D  is  not  possible  as  the 
depth  of  the  wake  centerline  from  the  free  surface  is  also  needed.  For  any  combination  of  (Re,Fr),  the 
linear  theory-based  evolution  curve  of  X^/D  could  lie  along  any  of  multiple  curves  whose  offset  depends 
on  the  depth  of  the  wake  core  (Figure  10).  In  more  practical  terms,  given  a  surface  observation  of 
wavelength,  the  passage  time  and  depth  of  the  wake-generated  body  cannot  be  inferred  simultaneously  ; 
at  least  one  of  these  quantities  must  be  known  a  priori.  One  must  know  the  depth  to  compute  how  much 
time  (in  Nt  units)  has  elapsed  since  the  passage  of  the  body.  Alternatively,  the  passage  time  must  be 
needed  to  estimate  the  depth  of  the  wave  source.  In  both  cases,  priori  knowledge  of  the  deeper  water 
stratification  is  also  necessary. 


10  10  10 


Figure  10:  Time  evolution  otXf/D  in  Nt  for  various  wake  depths,  Zsfi/D  E  [9,  IS,  36),  as 

predicted  by  our  linear  model  using  the  calibrated  virtual  origin  zO/D  =  0.62  for  waves  propagating 

at  0  =  45*^  to  the  vertical.  A  given  Xf/D  observed  at  the  surface  (marked  by  the 

horizontal  grey  line)  may  correspond  to  different  iV/ values  (where  the  vertical  grey 

arrows  point  at,  for  each  z^f/D  considered  in  the  figure)  of  the  submerged  turbulence 

which  may  be  located  at  various  depths. 

The  horizontal  wavelength  with  the  highest  probability  of  occurrence,  /  D ,  as  sampled  from  the  peak 
values  of  the  joint  p.d.f.’s  in  Figure  7,  is  shown  in  Figure  1 1,  For  a  given  Re,  an  approximate  scaling 
is  observed,  consistent  with  previous  numerical  and  experimental  studies.  In  terms  of  Reynolds  number 

dependence,  increasing  Re  by  a  factor  of  20  at  a  given  Fr  leads  to  a  50%  reduction  in  /  D .  A  first 
explanation  may  be  offered  in  terms  of  the  Ozmidov  scale  of  the  turbulence  at  the  onset  of  the  NEQ 
regime. 

Additional  work  has  been  performed  in  terms  of  the  potential  of  the  surface-reflecting  internal  waves  to 
generate  remotely  visible  slicks,  A  2-D  Fourier-based  advection  solver,  driven  by  the  horizontal 
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velocities  of  the  internal  waves  at  the  surface,  has  been  run,  initialized  by  an  initial  model  surfactant 
concentration  of /"q.  The  reflecting  wave  strain  fields  will  then  modulate  this  surfactant  distribution  and 
generate  regions  of  instantaneous  enhanced/reduced  surfactant  concentration,  corresponding  to  values 
ofr/ /"o  greater  or  larger  than  unity.  An  example  is  shown  in  Figure  12,  where  the  normalized  surfactant 
concentration  is  shown  along  with  the  normalized  horizontal  divergence  at  time  Nt=10  for  ^>=64  and 
/f^=10^  In  remote  sensing,  strain-driven  concentration  of/’  /  fo,  is  considered  to  be  remotely  visible 
slicks.  We  have  found  that  high  Re  wakes  have  a  pronounced  tendency  to  generate  such  slicks,  with  this 
tendency  increasing  with  Fr=64  as  shown  in  Figure  12b. 
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Figure  11:  Fr-  and  Re-dependence  of  the  most  energetic  wavelength  /  £)  observed  at 
the  sea  surface. 

Scaling  analysis  of  the  above  advection  equation  shows  that  the  normalized  surfactant  concentration 
perturbation  F'  /  Fq  should  scale  with  the  nonnalized  horizontal  divergence,  A  JN,  of  the  internal  wave 
field.  Figure  13  shows  that  the  LES  data  confirms  this  observation.  A  possible  leveling  off  of  the 
numerical  data  at  Re~\0^  and  Fr=64  requires  further  investigation  at  parameter  values  that  are 
inaccessible  within  current  computational  resources. 
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Figure  12:  Surface  snapshot  of  (a)  wave-induced  horizontal  divergence  z  and  (b)  straind riven 
surfactant  concentration  fields  at  Nt  =  70  for  and  F'r=64.  Concentration  of  F/  r'o> 

1,05  may  correspond  to  surface  slicks  that  can  be  detectable  via  remote  sensing.  Isolines 
with  r /Fo=  LOS  are  marked  in  grey,  delineating  possible  locations  of  slicks. 

Finally,  although  we  have  demonstrated  that  the  wake-radiated  internal  waves  do  not  undergo  nonlinear 
effects,  at  least  once  they've  left  the  near-field,  they  are  still  likely  to  be  subject  to  such  effects  upon  their 
reflection  at  the  free  surface.  This  hypothesis  is  driven  by  the  observation  of  nonlinear  effects  (generation 
of  harmonics  and  Lagrangian  mean  flow  formation)  we  have  observed  in  our  idealized  studies  of 
reflecting  2-D  internal  wave  beams.  To  this  end,  2-D  Lagrangian  particle  tracking  has  been  performed  on 
the  free  surface*  A  large  number  of  tracer  particles  are  inserted  onto  the  free  surface,  prior  to  the  arrival 
of  any  internal  wave  activity.  Particles  are  spaced  apart  by  O.IF.  By  virtue  of  the  free-slip  condition  (w=0 
at  the  impermeable  boundaty),  these  particles  cannot  be  move  downward  and  are  only  displace  on  the 
free  surface*  One  may  then  examine  the  statistics  of  lateral  particle  displacement  z(|y  |*  A  net  positive 
displacement  suggests  that  any  Lagrangian  mean  flow  due  to  reflecting  waves  is  sufficiently  powerful  to 
slow  push  particles  outward  away  from  the  wake  centerline. 
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Figure  13:  Peak  surfactant  perturbation  ratio  F  7  Fq  versus  the  peak  dimensionless  surface 
strain  |  z|/N  observed  in  all  six  simulations.  Only  F  >  0,  i,e,,  surfactant  enrichment, 
and  z  <  0,  i.e,,  flow  convergence,  are  considered.  The  characteristic  enrichment  ratio 
F  7  Fo=  0.05  (dash  line)  across  observable  slicks  as  commonly  proposed  in  oceanographic  field 
measurements  is  exceeded  by  all  cases  except  for  Fr=4  at  the  low  Reynolds  number,  A  scaling  of 
F  7  Fo  ^  /W (solid  line)  can  be  observed. 

Computational  restrictions  have  prevented  us  from  examining  the  Re-IQ^,  Fr=64  case*  However,  non- 
negligible  particle  drift  is  seen  for  Fr=4  and  16  at  this  higher  Reynolds  number  (Figure  14).  Prior  to  the 
arrival  of  the  most  energetic  waves  at  the  surface  (dashed  lines  in  Figure  14),  the  p,d*f.  of  lateral  particle 
displacements  centers  around  zero  with  some  symmetric  deviation*  After  the  most  intense  wave  impact 
(solid  lines  in  Figure  14),  a  biased  towards  non-negligible  particle  drifts  away  from  the  wake  centerline  is 
observed.  Given  that  nonlinear  effects  will  increase  with  both  Fr  and  i?e,  it  is  not  unlikely  that  this 
particle  can  be  as  large  as  0(D)  at  operationally  relevant  parameter  values*  For  an  underwater 
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submersible,  this  might  suggest  approximately  10m  displacements  of  surface  particulates  which  are 
clearly  visible  through  remote  sensing  technology. 


Our  findings  on  the  surface  manifestation  of  a  submerged  stratified  turbulent  wake  may  be  found  in  the 
article  by  Zhou  and  Diamessis  {J.  Fluid  Mech.  20 1 6). 


Figure  4:  PDF  of  lateral  particle  displacements  of  away  from  the  wake  centreline  at  (a) 

/?e=10*  and  Fr=4  and  (b)  /?e=I0*  and  Fr=16.  Dash  lines  denote  the  PDF  before  the  most  energetic 
waves  interact  with  the  surface,  and  the  solid  lines  the  PDF  after  the  most  energetic  wave 
impact. 


(Hi)  Very  high  Reynolds  number  stratified  wakes. 

Enabled  by  a  DoD-HPCMP  Frontier  grant  (PI:  Prof  Steve  de  Bruyn  Kops,  U.  Mass.,  Amherst).,  three 
simulations  at  /?e=4xl0^  and  Fr=4,  16  and  64  run  up  to  time  1,000  were  completed  on  DoD-HPC 
systems  using  a  dedicated  3,072  processor  partition.  These  three  simulations  are  at  the  highest  sphere- 
based  Reynolds  number  achieved  so  far.  Moreover,  they  effectively  provide  three  critically  important 
additional  points  in  the  (Re,Fr)  parameter  space  to  complement  existing  data  at  the  same  three  Fr  values 
and/?e=5xl0^and  10^ 

Each  7?e=4xl0^run  required  1024x1024x1189  grid  points  and  used  1024  cores  (512  MPl  processes). 

The  choice  of  resolution  was  made  to  enable  a  broader  dynamic  range  than  what  was  considered  in  the 
/?e=10^  case.  lOM  CPU  hrs  were  consumed  in  FY  2014  through  meticulous  simulation  supervision  by  the 
PhD  student  funded  by  this  grant.  An  approximately  20Tb  large  database  has  been  generated. 
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Figure  15:  Contour  plots  of  vertical  vorticity  (Uj  at  A'if=150  sampled  at  the  Oxy  horizontal  mid-plane 
for  simulations  at  Re  =  5>«10^,  10*  and  4x10*  respectively  (from  left  to  right)  and  Fr=4.  Sphere 
travels  from  left  to  right.  Distance  between  adjacent  ticks  on  the  axes  is  5/>. 

Figure  15  shows  example  stream-span  cuts  for  the  late  wake,  at  A^/=150.  The  late-time  pancake  vortices 
of  low  Re  stratified  wakes,  typical  of  the  laboratory,  are  replaced  quasi-horizontal  vortices  carrying 
significant  fine-structure  which  persists  at  times  as  late  at  /Vt=300  at  /ie=4x]0*.  This  turbulent  fine- 
structure  is  intimately  connected  to  secondary  Kelvin-Helmholtz  instabilities  visible  on  stream-depth 
cuts,  resulting  from  the  enhanced  buoyancy-induced  shear  at  higher  Reynolds  numbers. 

Following  the  pioneering  work  of  Spedding  (J.  Fluid  Meek  1997),  we  have  sought  to  identify  the 
transition  point  between  the  NEQ  and  Q2D  regimes,  as  manifested  by  a  change  in  the  decay  power  law  of 
the  mean  centerline  velocity  Uq.  at  this  point,  the  power  law  changes  from  a  -1/4  to  a  -3/4  exponent. 
Although  the  NEQ  regime  appears  to  last  consistently  longer  with  increasing  Re,  we  were  unable  to 
identify  a  distinct  transition  in  exponents  as  postulated  by  Spedding. 


Figure  16:  LEFT  PANEL:  Evolution  of  the  buoyancy  Reynolds  number,  Ret^RenFrH,  as  a  function 
of  time  in  buoyancy  units  for  ail  nine  stratified  wake  simulations.  Here  Re// and  Frff  are  local 
Reynolds  and  Froude  numbers  based  on  the  horizontal  integral  scale  and  characteristic  horizontal 
velocity  of  the  turbnlence.  Purple,  blue  and  red  line  colors  correspond  to  Re=4x  1 0* ,  1 0*  and  4x  ]  0* , 
respectively.  Solid,  dashed  and  dotted  lines  correspond  to  Fr=4,  1 6  and  64,  respectively.  RIGHT 
PANEL:  Buoyancy  Reynolds  number,  Ret,  rescaled  with  the  sphere-based  Reynolds  number,  Re. 
The  key  idea  is  that  the  time  where  Re/,  crosses  unity  scales  with  Re,  indicating  that  at  operationally 
relevant  Re,  the  NEQ  regime,  and  its  associated  stratified  turbulence,  can  persist  over  extremely 
long  times. 
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More  robust  information  on  the  transition  point  between  NEQ  and  Q2D  and,  therefore,  the  duration  of 
the  former  may  be  obtained  by  examining  the  buoyancy  Reynolds  number,  J?,  defined  in  terms  of 
characteristic  horizontal  velocities  and  the  horizontal  integral  lengthscale  of  the  stratified  turbulence, 
provides  further  insight  into  the  persistence  of  the  NEQ  regime  (Figure  16,  left  panel),  i?  is  a  measure  of 
the  separation  between  the  Ozmidov  (largest  scale  which  can  overturn  against  the  restoring  effect  of 
buoyancy)  and  the  Kolmogorov  scales  of  the  turbulence.  Once  it  drops  below  a  value  of  !,  stratified 
turbulence  is  suppressed  and  a  transition  to  the  02D  has  taken  place.  Consistently  with  figure  15,  this 
crossing  of  the  threshold  occurs  later  with  increasing  Re. 

Furthermore,  this  crossing  occurs  at  a  time  that  scales  with  Re  (Figure  16,  right  panel).  For  example,  if  at 
7?^”1 0^,  this  transition  occurs  at  Ni^]00^  at  the  operationally  relevant  value  of  the  transition  out 

of  the  NEQ  regime  would  take  place  at  0,000.  If  a  typical  value  of  oceanic  buoyancy  frequency 
ranges  between  0.01  and  0.1  rad/sec  several  questions  emerge  about  the  persistence  of  a  stratified 
turbulent  wake  in  the  field. 

Our  remaining  thrust  area  in  the  framework  of  high  Re  stratified  wakes  is  the  quantification  of  turbulent 
transport  coefficients,  i,e,,  horizontal  and  vertical  eddy  viscosities.  These  eddy  viscosities  are  integral  to 
rapid-turnaround  self-similarity  models  for  the  prediction  of  downstream  mean  wake  velocity  evolution. 
Vorticity  visualizations  (see  Figure  15)  indicate  that  the  horizontal  spreading  rates  of  higher  Re  stratified 
wakes  might  be  faster  than  their  low  Re  counterparts.  Moreover,  the  presence  of  persistent  stratified 
turbulence  extending  up  to  suggests  that  non-negligible  vertical  eddy  diffusivities  should  be 

present  well-after  the  critical  time  of  the  onset  of  the  NEQ  regime,  after  which  Meunier  et  al, 

{Phys.  Fluids  2006)  proposed  that  the  vertical  Reynolds  stresses  go  to  zero.  Such  a  proposal  is  clearly 
biased  by  a  reliance  on  low  Re  experimental  data. 

The  procedure  of  computation  of  the  horizontal  eddy  viscosity  is  illustrated  for  one  instant  {Nt=2)  for  the 
Fr=4  and  Re=lQ^  case  in  Figure  17*  Both  the  transverse  gradient  of  the  mean  velocity  and  the  horizontal 
Reynolds  stresses  are  averaged  along  the  homogeneous  x-direction.  For  values  of  horizontal  Reynolds 
stress  exceeding  a  threshold  value  (40%  of  the  max),  a  scatter  plot  of  the  above  two  quantities  is 
constructed,  A  linear  least  squares  fit  to  this  cloud  of  data  points  has  a  slope  which  is  equal  to  the 
horizontal  eddy  viscosity  of  the  stratified  wake  at  the  particular  time. 


Figure  17:  Left:  Streaiuwise-averaged  jz-trausect  of  the  lateral  gradient  of  the  mean  velocity  at 
for  the  and  Fr=4  (top)*  Equivalent  streamwise-averaged  transect  of  horizontal 

Reynolds  stress  field  Right:  Scatter  plot  of  these  quantities  for  Reynolds  stresses  that 

exceed  40%  of  the  maximum  value*  The  slope  of  the  linear  least  squares  fit  represents  a  measure  of 
the  eddy  viscosity,  R  is  the  correlation  coefficient. 
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up  to  Nt^2^  the  wake  behaves  effectively  likes  its  non-strati fied  counterpart*  Hence^  horizontal  and 
vertical  eddy  viscosities  are  equah  Beyond  this  point,  a  similar  procedure  to  the  above  is  used  to  compute 
the  vertical,  buoyancy*influenced,  eddy  viscosity*  This  procedure  is  the  object  of  ongoing  scrutiny  and 
discussion  (see  below). 

The  time  evolution  of  the  horizontal  eddy  viscosity,  for  all  9  cases,  is  shown  in  Figure  18.  Quite 
curiously,  the  eddy  viscosity  departs  from  the  typically  assumed  constant  value  at  Non- 
dimensionalizing  with  Uo^nd  Lh  ,  mean  profile  quantities,  does  not  seem  to  collapse  the  data  past 
i*e.  a  time  when  stratification  has  dominated  the  flow  (Figure  18a).  Moreover,  the  constant  eddy  viscosity 
values  proposed  in  independently  by  Bevilaqua  and  Lykoudis  (J*  Fluid Mech  1978)  and  the  classical 
textbook  of  Tennekes  and  Lumley  (1972)  do  not  seem  to  match  our  results.  The  former,  used  in  the 
stratified  wake  self-similarity  model  of  Meunier  et  aL  {Phys.  Fluids  2006),  shows  the  strongest 
disagreement.  Non-dimensional izing  with  the  square  root  of  the  maximum  horizontal  Reynolds  stress  and 
an  alternative  choice  of  characteristic  lengthscale  linked  to  the  shear  in  the  mean  profile,  proposed  by 
Tennekes  and  Lumley,  collapses  the  data  much  better  (Figure  18b)*  The  physical  implications  of  this 
improved  scaling  are  currently  under  investigation* 


Figure  18;  (a)  Horizontal  eddy  viscosity  for  all  9  wake  simulations,  scaled  with  the  mean  centerline 
velocity,  and  the  wake  half-width,  Lh*  (b)  Same  quantity  but  sealed  with  u*^  the  maximum 
horizontal  Reynolds  stress,  and  /*  an  effective  Reynolds  stress  profile  e- folding  scale  (see  the 
textbook  by  Tennekes  and  Lumley). 
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Figure  19:  Vertical  eddy  viscosity  for  all  nine  wake  simulations  scaled  with  mean  wake  centerline 
velocity  and  mean  profile  half-height.  The  correlation  of  the  least  squares  fit  (equivalent  to  what  is 
shown  in  Figure  14)  Is  lost  after  Vl^lO,  presumably  due  to  the  approach  used  in  computing  the 
vertical  eddy  viscosity  which  does  not  focus  on  physics  particular  to  stratified  turbulence. 

Results  for  the  vertical  eddy  viscosity,  using  mean-flow-quantity-based  scaling,  are  shown  in  Figure  19. 
The  particular  scaling  seems  to  collapse  the  results  from  all  9  cases  up  to  AftwlO.  Beyond  this  point,  any 
correlation  in  the  scatter  plots  equivalent  to  the  one  shown  in  Figure  1 7  is  lost  and  a  no  reliable  slope  of  a 
linear  squares  fit  can  be  computed  and,  therefore,  no  eddy  viscosity  estimate  is  possible. 

It  is  here  where  it  is  very  likely  that  the  vertical  eddy  viscosity  cannot  be  computed  in  a  manner 
analogous  to  its  horizontal  counterpart.  This  was  the  focal  point  of  the  Pi’s  discussions  with  Dr.  Patrice 
Meunier  during  his  ONR-funded  summer  visit  to  IRPHE  Marseille  in  France.  Vertical  transport  is 
dominated  by  r-iocalized  processes  (secondary  K-H  instabilities)  the  underlying  shear  of  which  is  not 
identical  to  that  of  the  mean  flow.  Hence,  one  must  examine  the  rms-value  of  the  vertical  shear  of  the 
horizontal  and  vertical  velocities  to  correctly  account  for  the  full  effect  of  buoyancy-driven  shear.  A 
conditional  sampling,  considering  local  values  of  the  Richardson  number  below  unity,  might  also  be 
necessary.  Finally,  it  is  likely  that  the  vertical  profile  might  not  necessarily  be  a  Gaussian  and  effectively 
the  solution  to  the  equations  of  high  Re  wake  mean  profile  evolution  may  not  be  self-similar.  In  other 
words,  the  modified  mean  profile  evolution  equations  may  only  be  solved  numerically.  These  issues  are 
the  focal  point  of  our  current  investigations. 

The  most  recent  update  on  our  work  on  high  Re  effects  in  stratified  wakes  may  be  found  in  the 
conference  article  by  Diamessis  and  Zhou  (2016). 


(iv)  Collaborative  DNS/LES  of  high  Reynolds  number  sti'atified  wakes  (Frontier- 
award  support) 

Beyond  generating  an  unprecedented  dataset  at  high  Re,  our  Frontier-project-driven  collaboration  with 
Prof,  de  Bruyn  Kops,  Prof.  Jim  Riley  (U.  Washington)  and  Dr.  Andreas  Muschinski  (NWRA)  has 
advanced  considerably  since  its  initiation  in  Fall  2013.  Apart  from  regular  weekly  discussions,  the  most 
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notable  accomplishment  has  been  Prof,  de  Bruyn  Kops  execution  of  a  subset  of  DNS  runs  initialized  by 
our  wake  data  at  Nt=3^  11,  21,41  and  6  at  Fr=4  and  ^e=l  0^  and  4x10^.  Each  run  is  conducted  for  a 
duration  of  Nl=2  which  allows  for  the  small-scale  portion  of  the  turbulence  spectrum,  unresolved  by  the 
LES,  to  be  “filled  up’^  Analysis  of  the  existing  DNS  results  by  a  postdoc  visiting  Prof  Riley  has  yielded 
significant  insight  on  the  non-trivial  process  of  identifying  the  wake’s  non-turbu lent/turbulent  interface  in 
a  stratified  environment,  the  spatial  distribution  of  turbulence  lengthscales  w/r  to  this  interfacial  zone  and 
the  presence  of  internal  wave-radiated  energy  that  is  comparable  to  the  turbulent  kinetic  energy 
dissipation  inside  the  wake.  This  work  was  published  in  Watanabe  et  al  (J.  Fluid Mech  2016). 


Ongoing/Future  Work 

Current  work  is  focused  on  completing  the  high  Re  scaling  analysis  of  our  stratified  wake  database  as 
presented  in  section  (iii)  of  this  report.  The  aim  is  to  provide  a  theoretical  foundation  for  the  observed 
scaling  and  a  first  predictive  parameterization  of  the  trajectory  of  a  stratified  turbulent  wake  in  the 
turbulent  Reynolds/Froude  number  database.  A  postdoctoral  researcher  funded  by  the  continuation  of 
this  grant,  arrived  at  Cornell  on  4/1/2016,  and  has  already  completed  systematic  work  on  quantifying  and 
parameterizing  the  interna!  wave  energy  fluxes  out  of  a  stratified  turbulent  wake.  In  parallel,  the  postdoc 
is  exploring  the  relevant  significance  of  these  fluxes  in  the  greater  wake  energy  budget  and  is  examining 
the  potential  for  sound  generation  at  high  Froude  numbers.  As  a  next  step,  the  postdoc  will  work  on 
optimizing  the  current  LES  code,  for  optimal  performance  on  Knight’s  Landing  Architectures  on  state- 
of-the-art  DoD  HPC  systems.  He  will  also  work  on  implementing  a  pilot-level  strategy,  currently 
implemented  in  MATLAB,  for  minimizing  spurious  divergence  at  spectral  subdomain  interfaces  in  the 
actual  Fortran  90/95  code.  This  strategy  will  allow  for  more  stable  wake  simulations  at  body-based 
Reynolds  numbers  as  high  as  2x  l  O^  which  have  never  been  attained  so  far. 


Personnel 

This  ONR  grant  has  supported  one  Ph.D.  student  for  the  entire  course  of  his  thesis.  Dr,  Qi  Zhou  received 
his  PhD  in  May  2015.  His  2-year  postdoc  at  Cambridge  U.  is  nearing  an  end.  In  September  2017,  he  will 
join  the  faculty  of  the  Department  of  Civil  and  Environmental  Engineering  at  the  University  of  Calgary 
in  Canada  as  an  Assistant  Professor.  The  postdoctoral  researcher  currently  supported  by  the  continuation 
of  this  grant  is  Dr.  Kristopher  Rowe,  who  was  granted  a  PhD  in  Applied  Mathematics  by  the  University 
of  Waterloo  in  March  2016. 
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